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Abstract. We study the phase behavior of solutions consisting of positive and negative ions of valence 
z to which a third ionic species of valence Z > z is added. Using a discretized Debye-Hiickel theory, 
we analyze the phase behavior of such systems for different values of the ratio Z = Z/z. We find, for 
Z > 1.934, a three-phase coexistence region and, for Z > 2, a closed (reentrant) coexistence loop at 
high temperatures. We characterize the behavior of these ternary ionic mixtures as function of charge 
asymmetry and temperature, and show the complete phase diagrams for the experimentally relevant cases 
of Z = 2 and Z = 3, corresponding to addition of divalent and trivalent ions to monovalent ionic fluids, 
respectively. 
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1 Introduction 

In many areas of Physics, Chemistry, and Biology elec- 
trostatic forces play an important role in determining sys- 
tem properties. This becomes particularly true when deal- 
ing with ionic fluids, i.e., fluids consisting of dissociated 
cations and anions, since in most cases the Coulomb inter- 
action is the dominant interaction between the particles. 
Ionic mixtures have been intensely debated not only be- 
cause of their interesting behavior (for recent reviews, see 
||l|,||||,§,||), but also because a better understanding of 
such systems is the first step towards a more rigorous theo- 
retical treatment of more complicated systems like charged 
. colloids and polyelectrolytes. 

It has been shown experimentally that ionic solutions 
exhibit liquid-liquid phase coexistence, terminated by a 
' critical point, similar to ordinary (nonionic) fluids B. The- 
[ oretically, such a demixing transition has been rationalized 
. in terms of an effective attraction between the ions|0,|| 
' which is due to charge screening, first introduced by Debye 
and Hiickeip]. For simple inorganic salts, such as NaCl, 
this critical point occurs at temperatures above 3000 Kg 
and thus precludes precise measurements. More detailed 
experiments became available with the advent of large or- 
ganic ions, which show critical points at temperatures of 
414 K|i^ and, more recently, at 317 K0. 

The aforementioned examples correspond to binary 
mixtures of cations and anions of the same valence. In 
this article we will look at three-component ionic mix- 
tures, i.e., we will consider solutions of z:z salts to which 
a third component with valence Z is added (keeping to- 
tal charge neutrality). This situation often occurs in real 
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systems, deserving for this reason a closer look. For in- 
stance, in the cellular medium several ionic species with 
different valences (like Cl~, Na"*" and Ca^+) are simul- 
taneously present. As another example, the presence of 
charged impurities can change dramatically the behavior 
of charge-stabilized colloidal suspensions |13] . With this in 
mind, adding a third component to an ionic fluid is the 
natural extension of previous work and, as we will show, 
the resulting phase behavior is quite rich: it is our hope 
that this will motivate experiments testing our predic- 
tions. 

The theoretical tools that are usually employed to treat 
ionic fluids are the Debye-Hiickel theory^ (or some modi- 
fied version of it) and liquid state theory with some special 
approximation ||l^ like the hypernetted chain (HNC) or 
the mean-spherical approximation (MSA). For instance, 
Caccamop^ and Caccamo and Malescio p5[| studied ternary 
mixtures with Z = Z/z = 2, the former using MSA (for 
low-density mixtures) , while the latter made a comparison 
between MSA, HNC and Monte Carlo results (at high salt 
densities). Also, three-component mixtures with Z = 2 
and 3 have been studied by Kenkare et al. |lq| , where MSA 
was used to compute the Coulombic contribution to the 
free energy. Their phase diagrams, although presented in 
a constant-pressure ensemble (in contrast to our phase di- 
agrams, where we do not constrain the osmotic pressure 
to constant values) , seem to be generally in accordance 
with our results |T7|,[T^. 

In general, the comparison between computer simula- 
tions and liquid state theory shows that the latter leads 
to quite accurate predictions for the thermodynamic vari- 
ables of the system. The same cannot be said about the 
Debye-Hiickel theory (cf. Table 1): although it should be 



Andre G. Moreira and Roland R. Netz: Phase behavior of three-component ionic fluids 



exact in the limit of vanishing densities, it quickly deviates 
from the results obtained by computer simulations as one 
increases the salt density. However, if one is interested in 
low-density systems, Debye-Hiickel has the advantage of 
offering a very simple free energy that, although not very 
accurate, can still be used to study the phase behavior 
of ionic systems without considerable numerical work (in 
contrast to most approximations used in liquid state the- 
ory). Besides, Debye-Hiickel has been shown ||l{|po| to be 
a Gaussian theory that can be improved, with the help of 
field-theoretic methods, through a systematic expansion of 
the free energy in cumulants of a fluctuating field (which 
can be related to the electrostatic potential). For these 
reasons, we will use here the free energy derived from the 
Debye-Hiickel theory to describe our system. 

Our article is presented as follows: in section 2 we 
briefly introduce the theoretical background of our work, 
in section 3 we present the resulting phase diagrams for 
various three-component systems and fully characterize 
their phase behavior as a function of Z and temperature 
and finally section 4 contains some concluding remarks. 



2 Theoretical background 

Let us assume a symmetric ionic system with cations and 
anions of valence z to which a third ionic species of valence 
Z is added. Let us also assume that the ions interact only 
via the Coulomb force, or that any solvophobic interaction 
between the ions are weak when compared to the Coulomb 
interaction. In this case, following what has been derived 
in Ref. |jl9, the Debye-Hiickel contribution to the bulk 
free energy per volume is 
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where the electrostatic self-energy is subtracted and the 
momentum integral goes over a cube of length lix j a. This 
accounts for the lattice (with spacing a) , implementing in 
an approximate way the ionic hard cores. The screening 
length K~^ is defined by 



47r^i3(z2c++z2c_ +Z2c), 



(2) 



where the concentrations of the positive and negative ions 
of valence z are denoted by c+ and c_ , and the concentra- 
tion of the third component (with valence Z) is denoted by 
c. The Bjerrum length Ib = e^ /{A:Tt eksT) (in SI units) 
is the distance at which the electrostatic energy between 
two elementary charges equals the thermal energy ks T. 

Since the integrand in Eq. (|l|) is isotropic, we distort 
the integration volume to a sphere and obtain after a 
straightforward integration 
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In the limit a ^ this reduces to /dh — —k.^/12tt, 
which is the well-known Debye-Hiickel limiting law. The 



full (Helmholtz) free energy density contains also the ideal 
entropy of mixing and reads 
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Global charge neutrality implies that Zc + zc^ = zc_ , 
where, without loss of generality, the third component of 
valence Z is assumed to be positively charged. This allows 
the elimination of one of the concentrations from the free 
energy in Eq. (jj) : in the following, we choose to eliminate 
c_. For representing our results, we introduce the usual 
rescaled temperature T* = a/lsz'^ and the volume frac- 
tions 0+ = a^c-\- and (p = a^c (notice that the usual defi- 
nition of volume fraction — normally represented as 77 — is 
related to our definition through = 6 ri/w) . 

If one of the three species is absent, the system reduces 
to a two-component system, which has been studied (us- 
ing the same free energy) in Ref. |13|. If </> = 0, we obtain 
the lower limiting critical temperature at T^i^^ — 0.1776 
with a critical density of positive ions of ^llq — 0.0209; 
if 0+ = 0, we obtain the upper limiting critical temper- 
ature at T*^c_= 0.1776 Z and 0ulc = 0.0418/(1 + Z). 
Clearly, since Z > 1, we have T^i^^. > T^^fj. These tem- 
peratures will be used as "guidelines" in the presentation 
of the phase behavior of three-component systems. 

In order to obtain the phase diagrams, we calculate 
the chemical potentials fi and /i"*" from the free energy 
Eq. (H) by taking derivatives with respect to c and c+ , re- 
spectively. We also calculate the osmotic pressure p, which 
can be easily shown to be p = —f + c^ + c+^j,'^ (in units of 
ksT). We look for instabilities in the free energy through 
the sign of the determinant of the (symmetric) matrix 



(5) 



in the plane (c, c+) (or equivalently in the plane (0, 0^)) at 
fixed values of T* . The global convexity of the free energy 
implies that det J should be always positive: when this 
is not verified, one is in a region where phase separation 
occurs. The coexistin g p hases are then found through the 
Maxwell construction |2 1| , i.e., one has to find the c and c+ 
in each of the coexisting phases that lead to values of ji, 
fi~^ and p that are same for all phases in coexistence. This 
procedure yields the phase diagrams that we now present. 



3 Phase diagrams 

In Fig. 1 we present the phase behavior for Z — 1.8, rep- 
resentative of small charge ratios Z < 1.934. For T* > 
"^ULC ~ 0.3196 no phase separation occurs. In Fig. la we 
show, to the left, the phase diagram in fugacities (which 
are defined as the exponential of the chemical potentials) 
and, to the right, the phase diagram in volume fractions 
for T* = 0.2500, representative of the range T£lc < T* < 
-^ULC- There is one critical point at the terminus of a phase 
coexistence line which emanates from the axis defined by 
0+ = (or e^^ = 0). At T* = T^^f. this critical point hits 
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Fig. 1. Phase diagrams at constant reduced temperature 
T* = a/lsz^ for a three-component system with a valence ratio 
Z = 1.8 — representative of systems with Z = Z/z < 1.934 — as 
function of e'' and e^ (the fugacities of .Z-valent ions and of 
positive z-valent ions, respectively) and as function oi (j> = a^ c 
and (j)^ = a'^ c+ (the volume fractions of Z-valent ions and 
of positive z-valent ions). Filled circles denote critical points. 
Plots: (a) T* = a/dsz^ ~ 0.2500 (representative of the tem- 
perature range T^^c < T* < T{,^c)\ (b) T* = 0.1754 (repre- 
sentative for T* < T{^i^q). 



the axis = (or e'^ = 0) at (?!>llc, and for lower temper- 
atures there is single coexistence line connecting the axes 
(/>"*" = and (j) = 0. This is demonstrated in Fig. lb, where 
we show a phase diagram for T* = 0.1754. 

The behavior is dramatically altered for larger charge 
asymmetry, viz. Z > 1.934. One important difference in 
the phase topology is that such systems always have two 
critical end-points, one at temperature T^cep (upper crit- 
ical end-point) and another one at T^c^p (lower critical 
end-point). In the temperature range between these two 
critical end-points the systems exhibit three phase coex- 
istence. 

In Fig. 2 we present the phase behavior for Z = 2 
corresponding, for instance, to a mixture of a divalent 
cations with monovalent co- and counter-ions. For high 
temperatures, T* > T{Jj^q = 0.3552 no phase separation 
occurs. For a temperature range T^Q^p < T* < Tulq 
there is a single critical point, as demonstrated in Fig. 2a 
for T* = 0.2070; the phase diagram is similar to the one 
shown in Fig. la. The upper critical end-point is located 
at ^ucEP = 0.1808, below which a second critical point 
exists. The resulting phase behavior for T^Q^p < T* < 
^UCEP is similar to Fig. 2b, which depicts the phase di- 
agram for T* — 0.1800, and where we see two phase co- 
existence lines, each terminating at a critical point and 
merging at a triple point. The region around this triple 
point is highlighted at the phase diagram in volume frac- 
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Fig. 2. Phase diagrams at constant T* = u/Ibz^ for 2 = 2. 
The notation is the same as in Fig. 1. Shaded areas in the 
volume fraction representation denote three phase coexistence 
zones. Plots: (a) T* = 0.2070 (T^cep < T* < T^j^c); (b) 
T* = 0.1800 (Tlcbp < T* < TucEp) (only shown close to the 
three phase coexistence zone); (c) T* = 0.1786 (Tllc < T* < 
TlcEp); (d) T* = 0.1429 (T* < T^lc)- 



tions (Fig. 2b, right). At T* = T^^ep = 0.1796 one of 



the critical lines present for TI^q^p < T* < Tucep termi- 
nates at the lower critical end-point. For T^lj-, < T* < 
-^LCEP ^e are thus left again with a single critical point; 
a representative phase diagram is shown in Fig. 2c for 
T* = 0.1786. Finally, for T* < T^^^ = 0.1776, no critical 
point is observed; Fig. 2d shows such a phase diagram for 
T* = 0.1429, which is similar to the one in Fig. lb. 

The complete phase behavior is summarized in Fig. 3, 
where we plot the various critical temperatures as a func- 
tion of the valence ratio Z. In Fig. 3a we show, in the 
vicinity oi Z = 2, the upper critical end-point temperature 
'^UCEP (broken line), the lower critical end-point temper- 
ature Tlcep (dotted- broken line), and the lower limiting 
critical temperature T^^i^^. (solid line). The two hnes of 
critical end-points meet at a double critical end-point at 
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Fig. 3. (a) Plot of the upper critical end-point temperature 
TucEP (broken line), the lower critical end-point temperature 
Tlcbp (dotted- broken line), and the lower limiting critical tem- 
perature Tllc (solid line) as a function of the valence ratio Z 
in the vicinity of the double critical end-point (DCEP). (b) 
Plot of the lower and upper limiting critical temperature, the 
maximal critical temperature (Tcmax) smd the lower critical 
end-point temperature as a function of Z\ notice the widen- 
ing of the three phase zone and the growth of Tcmax ^s Z 
increases. 



Z = 1.934, below which no triple point occurs. On a larger 
scale (Fig. 3b) it becomes clear that the triple phase tem- 
perature range widens as Z increases (note that lu^gp 
and the double critical end-point cannot be distinguished 
from T£l(-, at this scale). We also plot the maximal tem- 



perature at which a critical point occurs, T^ 



ioT Z > 2 satisfies T^max > ^ulc 



CMAX' 



which 



This leads to closed 
coexistence loops, as is demonstrated in the following. 

In Fig. 4 we show the phase behavior for Z = 3. 
In Fig. 4a we plot the phase diagram for T* = 0.5747, 
which is below the maximal critical temperature Tqjvjax ~ 
0.5771 and above the upper limiting critical temperature 
0.5328. Notice the presence of two critical points 



ULC 



as depicted in Fig. 4b for T* = 0.5155. In this system the 
upper critical end-point occurs at a temperature very close 
to the lower limit critical temperature (this is already clear 
in Fig. 3); more precisely, T]*lc = 0.1775656 and T, 



UCEP 



In Fig. 4c, where T* = T£lc+0.7x 10 



T£lc + 1.0x10- 
(^LLC < T* < T^cEp)^ ^^ observe a triple point and two 
critical points. Since the lower critical end-point occurs at 
a lower temperature than the lower limiting critical point 
(cf. Fig. 3), we obtain in the range Tl^^p < T* < Tl^^. 
one critical point and a triple point, as shown in Fig. 4d for 
T* = 0.1667. Finally, for T* < T^^^p = 0.1163, all crit- 
ical points have disappeared and a coexistence line con- 
nects the two phase diagram axes, as shown in Fig. 4e for 
T* = 0.1000. 



bounding a closed coexistence loop. For T^Qp^p < T* < 
-^ULC ^^ ^^^ ^ phase topology with a single critical point. 



4 Conclusion 

For small valence ratios {Z = Z/z < 1.934) we have shown 
that the the two critical points obtained for the limiting 
cases (i.e., one of the two positively charged species is ab- 
sent) are joined by a continuous critical line. For larger 
valence ratios {Z > 1.934) we have obtained more compli- 
cated phase diagrams, with a line of triple points bounded 
by critical end-points and closed coexistence rings at T > 
^ULC ior Z > 2. We have also made a global analysis of 
three-component ionic mixtures, where their phase behav- 
ior is characterized as a function of charge asymmetry (up 
to ^ ~ 5) and temperature. 

Denoting by T^^j^q the lower limiting critical tempera- 
ture of the symmetric ionic system, we have shown that 
this triple line extends, for Z = 2, over a temperature 
window of 1.012 T^c < T* < 1.018 T]*lc and, for Z = 3, 
over 0.66T]*Lc < T* < Tl^^. This means that for a 1 : 1 
organic salt with a critical point at T ~ 400 K this tem- 
perature window would have a 2.4 K width when adding 
divalent ions and a 136 K width when adding trivalent 
ions. The triple line should in principle be observable with 
a suitable choice of mono-, di- and trivalent organic ions. 
However, one should be careful when examining the num- 
bers presented here. Although the Debye-Hiickel theory 
successfully predicts for ionic systems the existence of a 
liquid- liquid critical point and the high asymmetry of their 
phase diagrams, it has some limitations, as we now discuss. 

One problem is the fact that the theory is quite sensi- 
tive to the way that the hard cores are treated: we chose 
to include them through a lattice of spacing a. Alterna- 
tively, one could add to the free energy Eq. (^ an extra 
term to account for the hard-core repulsion (for instance a 
free- volume approximation term) and get the electrostatic 
contribution /dh either by integrating Eq. (nl) without the 
cut-off (which results immediately in the Debye-Hiickel 
limiting law) or by doing a "charging process" after solv- 
ing of the linearized Poisson-Boltzmann equation with the 
appropriate boundary conditions for the hard cores [y,|2j, 
p3| . These methods yield similar phase diagrams, but with 
different values of temperature and density for the critical 
point. The use of a lattice is the most convenient for the 
field-theoretic treatment, although it is certainly not the 
best wayP4|- This, as well as some simplifications done 
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Fig. 4. Phase diagrams at constant T* = a/lBZ^ for Z — Z. 
The notation is the same as in Figs. 1 and 2. Plots: (a) T* = 
0.5747 (T^Lc < r* < T5max); (b) T* = 0.5155 (T^cep < 
T* < T^Lc); (c) T* = 0.1775663 (Tl^c < T* < T^cep) (in 
volume fraction the three phase coexistence is represented by 
a dashed line); (d) T* = 0.1667 (T^cep < T* < Tllc); (e) 
T* = 0.1000 (T* < Tl*cep). 



during the calculations (like performing the integral in 
Eq. (n^) in the sphere instead of the initial cube of size 
l-K ja) are expected to affect the final results. 

The direct comparison of the internal energy per par- 
ticle that follows from Eq. (Eh, viz. 
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Table 1. Comparison between (mirms) the excess internal en- 
ergy per particle from Monte Carlo simulations (MC), hyper- 
netted chain (HNC), mean spherical approximation (MSA), 
Debye-Hiickel limiting law (DH) and Eq. (H) (<^+ = a^ c+, 
as defined in Section 2). In the flrst part of the table ("bi- 
nary") the comparison is done for a simple 1 : 1 electrolyte 
with T* = 0.59, where the data for MC, HNC and MSA are 
taken from Table H of Ref. [g5| . In the second part of the table 
("ternary"), the same comparison is done for a three compo- 
nent mixture with Z = 2, T* — 0.294 and (f> — 0"|^ the data 
for HNC and MSA are taken from Table I of Ref. Q. 



with results from other approaches is shown in Table 1 (in 
Eq. (|6|), K is the inverse screening length defined in Eq. (y), 
a is the lattice spacing, corresponding to the ionic diam- 
eter, and (/)+ is the volume fraction of the positive z = 1 
ions). In the first part of this table ("binary"), we com- 
pare the values of u from Monte Carlo simulations, HNC, 
MSA [HI and Debye-Hiickel limiting law (DH) for a sim- 
ple 1 : 1 ionic mixture (where (j) = 0, i.e., the Z-valent ions 
are not present) . Although our results are better than the 
pure DH results, it is clear that our theory needs improve- 
ment at higher volume fractions. In the second part of ta- 
ble 1 ("ternary"), the same comparison is done for a three- 
component mixture with 2 = 1 with results obtained by 
CaccamoflJI (notice that the value of </!)+ is in the region 
where our phase diagrams are shown). This shows that our 
theory underestimates the value of the internal energy of 
the system, which indicates that, beyond the problem with 
the hard cores, a more refined free energy is needed in or- 
der to capture correctly all features of the ionic fluids [Em. 
This idea is reinforced by the discrepancy between the 
critical parameters for a 1 : 1 salt obtained through MC 
simulations [|3 (which lie in the range T^ = 0.049-0.070 
and 0j = 0.013-0.040) with the ones obtained with the 



free energy used here (T^ 



0.178 and (/.+ 



0.021). 



With this in mind, one should regard the results we 
obtained as being of "zeroth-order." As previously men- 
tioned, the free energy used here is the lowest-order result 
following a systematic expansion. Changes in the phase 
behavior are expected as one takes into account higher 
order terms in the expansion, but as shown in Ref. [ |l9| 
for Z : z mixtures, the inclusion of terms up to the eighth 
cumulant in the fluctuating field affects the phase dia- 
grams appreciably only for values oi Z = Z/z > 5. For 
this reason, and in spite of the limitations of the Debye- 
Hiickel theory, we assume that the topology of the phase 
diagrams shown here are correct. 
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